%cd ' c:\data\restat\programs\matlab';

load 'c:\data\restat_12_9\analysis_data\syruplogitdata';


%set seed
randn('state',32);
rand('state',32);
%data=st_dataset;


id=data(:,1);
s_jt=data(:,2);
%s_jt=s_jt*8;
x1=data(:,3:20);
%nx1=[data(:,2) data(:,4:17)];
clear data;
load 'c:\data\restat_12_9\analysis_data\syrupblpiv';
iv=data(:,2:10);

brand=floor(id/1000000);
company=floor(brand/1000);

pre_merger_dum=cr_dum(company);
own_dum=pre_merger_dum(1:9,:);



%load 'u:\user11\mw562_rs\My Documents\restat\analysis_data\outside4';
load 'c:\data\restat_12_9\analysis_data\outside4';
mean_price=data(:,1);
mean_share=data(:,2);
meany=log(reshape(mean_share,9,47))-log(ones(9,47)-kron(ones(9,1),sum(reshape(mean_share,9,47))));


IV=[iv x1(:,2:18)];
clear iv data;

clear data;

          
horz=['OLS    IV'];
vert=['constant  ';
      'price     ';
      'private   ';];




nmkt =47*9;  
nbrn =9;     
cdid=kron([1:nmkt]',ones(nbrn,1));
cdindex = [nbrn:nbrn:nbrn*nmkt]';       


% create weight matrix
invA = inv([IV'*IV]);


% compute the outside good market share by market
temp = cumsum(s_jt);
sum1 = temp(cdindex,:);
sum1(2:size(sum1,1),:) = diff(sum1);
outshr = 1.0 - sum1(cdid,:);
outshr=outshr.*(outshr>=0);
y = log(s_jt) - log(outshr);

%grand means for elasticity calculations


mean_price_atavg=mean(reshape(x1(:,1),9,9*47)')';

mean_share_atavg=mean(reshape(s_jt,9,9*47)')';

%ols logit
beta_ols = inv(x1'*x1)*x1'*y;
res = y - x1*beta_ols;
%vcov_ols = 1/size(x1,1)*res'*res*inv(x1'*x1);
vcov_ols = inv(x1'*x1)*(x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1))*inv(x1'*x1);
se_vcov_ols=sqrt(diag(vcov_ols));
%NEWEY WEST SE'S HERE%

%Need to get lagged value of xe

x1_lag1=x1(1:3807-423,:);
x1_lag2=x1(1:3807-423*2,:);
x1_lag3=x1(1:3807-423*3,:);


x1_temp1=x1(423+1:3807,:);
x1_temp2=x1(423*2+1:3807,:);
x1_temp3=x1(423*3+1:3807,:);


res_lag1=res(1:3807-423,:);
res_lag2=res(1:3807-423*2,:);
res_lag3=res(1:3807-423*3,:);

res_temp1=res(423+1:3807,:);
res_temp2=res(423*2+1:3807,:);
res_temp3=res(423*3+1:3807,:);


G_0=x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1);
G_1=x1_lag1'.*(ones(size(x1_lag1,2),1)*res_lag1')*(res_temp1*ones(size(x1_temp1,2),1)'.*x1_temp1);
%Need to get twice lagged value of xe
G_2=x1_lag2'.*(ones(size(x1_lag2,2),1)*res_lag2')*(res_temp2*ones(size(x1_temp2,2),1)'.*x1_temp2);
%Need to get three times lagged value of xe
G_3=x1_lag3'.*(ones(size(x1_lag3,2),1)*res_lag3')*(res_temp3*ones(size(x1_temp3,2),1)'.*x1_temp3);



vcov_ols_newey = inv(x1'*x1)*(G_0+.5*(G_1+G_1'))*inv(x1'*x1);
se_vcov_ols_newey=sqrt(diag(vcov_ols_newey));
clear x1_temp1 x1_temp2 x1_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 G_0 G_1 G_2 G_3;



% IV regressions
beta_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*IV'*y;
res = y - x1*beta_IV;
IVres = IV.*(res*ones(1,size(IV,2)));
b = IVres'*IVres;
vcov_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*IVres'*IVres*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV=sqrt(diag(vcov_IV));


%Newey West S.E.'s here

IV_lag1=IV(1:3807-423,:);
IV_lag2=IV(1:3807-423*2,:);
IV_lag3=IV(1:3807-423*3,:);


IV_temp1=IV(423+1:3807,:);
IV_temp2=IV(423*2+1:3807,:);
IV_temp3=IV(423*3+1:3807,:);


res_lag1=res(1:3807-423,:);
res_lag2=res(1:3807-423*2,:);
res_lag3=res(1:3807-423*3,:);

res_temp1=res(423+1:3807,:);
res_temp2=res(423*2+1:3807,:);
res_temp3=res(423*3+1:3807,:);

IVres_lag1 = IV.*(res*ones(1,size(IV,2)));
IVres_lag2 = IV.*(res*ones(1,size(IV,2)));
IVres_lag3=IV_lag3.*(res_lag3*ones(1,size(IV_lag3,2)));

G_0=IVres'*IVres;
G_1=IV_lag1'.*(ones(size(IV_lag1,2),1)*res_lag1')*(res_temp1*ones(size(IV_temp1,2),1)'.*IV_temp1);
G_2=IV_lag2'.*(ones(size(IV_lag2,2),1)*res_lag2')*(res_temp2*ones(size(IV_temp2,2),1)'.*IV_temp2);
G_3=IV_lag3'.*(ones(size(IV_lag3,2),1)*res_lag3')*(res_temp3*ones(size(IV_temp3,2),1)'.*IV_temp3);
vcov_IV_neweywest = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*(G_0+.5*(G_1+G_1'))*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV_neweywest=sqrt(diag(vcov_IV_neweywest));
clear IV_lag1 IV_lag2 IV_lag3 IV_temp1 IV_temp2 IV_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 IVres_lag1 IVres_lag2 IVres_lag3;




%Approximate first-stage F-stat, leave out vcov d.f. adjustment as N is large relative to number of regresors...
b_1st=inv(IV'*IV)*IV'*x1(:,1);
res=x1(:,1)-IV*b_1st;
vcov_1st=inv(IV'*IV)*(IV'.*(ones(size(IV,2),1)*res')*(res*ones(size(IV,2),1)'.*IV))*inv(IV'*IV);
se_vcov_1st=sqrt(diag(vcov_1st));
R=[eye(9) zeros(9,17)];
q=9;
F=(R*b_1st)'*inv(R*vcov_1st*R')*(R*b_1st)/q;

alpha=beta_IV(1);
pre_merger_dum=cr_dum(company);
for t=1:47
        shares=mean_share((t-1)*9+1:t*9);
 	    mp=mean_price((t-1)*9+1:t*9);
        own_dum=pre_merger_dum((t-1)*9+1:t*9,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*9+1:t*9,:)=o2-o1;
        elas((t-1)*9+1:t*9,:)=omega((t-1)*9+1:t*9,:).*(1./shares*mp');
        markup((t-1)*9+1:t*9,1)=-inv(omega((t-1)*9+1:t*9,:).*(own_dum*own_dum'))*shares;
end
ivelas=median(reshape(elas',9,9,47),3)'

mc_hat=mean_price-markup;


company_post=company.*(1-(company==4))+(company==4)*3;


post_dum=cr_dum(company_post);
%post_dum=cr_dum(company);

options=optimset('Tolfun',1e-12,'TolX',1e-12,'Display','on');
ivflag=zeros(47,1)
first=1;


for i=1:47
	last=i*9;
    expall=exp(meany(:,i)-alpha*mean_price(first:last));
	[p_star(first:last,1),y,ivflag(i,1)]=fsolve(@logitbert_eq,mean_price(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end

ivc=reshape(100*((p_star-mean_price)./mean_price),9,47)';
ivchange=median(ivc((ivflag==1),:))

mcflag=zeros(47,1);
first=1;
for i=1:47
	last=i*9;
    %will redefine shares so they're on demand curve given post_p within
    %mc_logitbert
    expall=exp(meany(:,i)-alpha*mean_price(first:last));
    post_p=mean_price(first:last).*[1+0.0052;1+.0118;1+0.0125;1+0.0156;1+0.0272;1+0.0273;1-0.0108;1+0.003;1+0.0108];
    [mc_post(first:last,1),y,mcflag(i,1)]=fsolve(@mc_logitbert,mc_hat(first:last),options,expall,alpha,... 
    post_p, post_dum(first:last,:));
	first=last+1;
end

ivmc_change=reshape(100*((mc_post-mc_hat)./mc_hat),9,47)';
ivmcchange=median(ivmc_change((mcflag==1),:))




alpha=beta_ols(1);
pre_merger_dum=cr_dum(company);
for t=1:47
        shares=mean_share((t-1)*9+1:t*9);
 	    mp=mean_price((t-1)*9+1:t*9);
        own_dum=pre_merger_dum((t-1)*9+1:t*9,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*9+1:t*9,:)=o2-o1;
        elas((t-1)*9+1:t*9,:)=omega((t-1)*9+1:t*9,:).*(1./shares*mp');
        markup((t-1)*9+1:t*9,1)=-inv(omega((t-1)*9+1:t*9,:).*(own_dum*own_dum'))*shares;
end

mc_hat=mean_price-markup;


company_post=company.*(1-(company==4))+(company==4)*3;


post_dum=cr_dum(company_post);
%post_dum=cr_dum(company);

options=optimset('Tolfun',1e-5,'TolX',1e-20,'Display','on');
olsflag=zeros(47,1);
first=1;
for i=1:47
	last=i*9;
    expall=exp(meany(:,i)-alpha*mean_price(first:last));
	[p_star(first:last,1),y,olsflag(i,1)]=fsolve(@logitbert_eq,mean_price(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end
olselas=median(reshape(elas',9,9,47),3)'
olsc=reshape(100*((p_star-mean_price)./mean_price),9,47)';
olschange=median(olsc((olsflag==1),:));



mcflag=zeros(47,1);
first=1;
for i=1:47
	last=i*9;
    %will redefine shares so they're on demand curve given post_p within
    %mc_logitbert
    expall=exp(meany(:,i)-alpha*mean_price(first:last));
    post_p=mean_price(first:last).*[1+0.0052;1+.0118;1+0.0125;1+0.0156;1+0.0272;1+0.0273;1-0.0108;1+0.003;1+0.0108];
    [mc_post(first:last,1),y,mcflag(i,1)]=fsolve(@mc_logitbert,mc_hat(first:last),options,expall,alpha,... 
    post_p, post_dum(first:last,:));
	first=last+1;
end

mc_change=reshape(100*((mc_post-mc_hat)./mc_hat),9,47)';
olsmcchange=median(mc_change((mcflag==1),:))















diary output
'Table 11, Appendix' 
olselas
'Table 12, Appendix' 
ivelas
'Table 2, Column 6'
olschange=median(olsc((sum(reshape(mean_share,9,47))<1)',:))
'Table 2, Column 7'
ivchange=median(ivc((sum(reshape(mean_share,9,47))<1)',:))
olsmcchange=median(mc_change((sum(reshape(mean_share,9,47))<1)',:))
ivmcchange=median(ivmc_change((sum(reshape(mean_share,9,47))<1)',:))

diary off



%%NEED TO CHECK VCOV MATRIX!!!%

%number of bootstrap draws
sims=1000;
olspchange=zeros(sims,9);
olsmelas=zeros(9,9,sims);
olsflag=zeros(sims,47);
for i=1:sims
    bs_beta_ols=beta_ols+chol(vcov_ols_newey)'*randn(size(beta_ols,1),1);
    alpha=bs_beta_ols(1);
    markup_multi = zeros(9*47,1);
    elas = zeros(9*47,9);
    omega = zeros(9*47,9);

    % compute demand elasticities


    o1=alpha*mean_share_atavg*mean_share_atavg';
    o2=diag(alpha*mean_share_atavg);
    omega=(o2-o1);
    
for t=1:47
        shares=mean_share((t-1)*9+1:t*9);
     	mp=mean_price((t-1)*9+1:t*9);
        own_dum=pre_merger_dum((t-1)*9+1:t*9,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*9+1:t*9,:)=o2-o1;
        elas((t-1)*9+1:t*9,:)=omega((t-1)*9+1:t*9,:).*(1./shares*mp');
        markup((t-1)*9+1:t*9,1)=-inv(omega((t-1)*9+1:t*9,:).*(own_dum*own_dum'))*shares;
end
olsmelas(:,:,i)=median(reshape(elas',9,9,47),3)';
mc_hat=mean_price-markup;
company_post=company.*(1-(company==4))+(company==4)*3;
post_dum=cr_dum(company_post);

 options=optimset('Tolfun',1e-5,'TolX',1e-10,'Display','on');

 first=1;
 flag=zeros(1,47);
for j=1:47
	last=j*9;
    expall=exp(meany(:,j)-alpha*mean_price(first:last));
	[p_star(first:last,1),y,x]=fsolve(@logitbert_eq,mean_price(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
    flag(:,j)=x==1;
	first=last+1;
end

olsflag(i,:)=flag;
c=reshape(100*((p_star-mean_price)./mean_price),9,47)';
check=(olsflag(i,:)'+(sum(reshape(mean_share,9,47))<1)'==2);
olspchange(i,:)=median(c(check,:));

end



%number of bootstrap draws
sims=1000;
ivpchange=zeros(sims,9);
ivmelas=zeros(9,9,sims);
ivflag=zeros(sims,47);
for i=1:sims
    bs_beta_IV=beta_IV+chol(vcov_IV_neweywest)'*randn(size(beta_IV,1),1);
    alpha=bs_beta_IV(1)
    markup_multi = zeros(9*47,1);
    elas = zeros(9*47,9);
    omega = zeros(9*47,9);

    % compute demand elasticities


    o1=alpha*mean_share_atavg*mean_share_atavg';
    o2=diag(alpha*mean_share_atavg);
    omega=(o2-o1);
    
for t=1:47
        shares=mean_share((t-1)*9+1:t*9);
     	mp=mean_price((t-1)*9+1:t*9);
        own_dum=pre_merger_dum((t-1)*9+1:t*9,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*9+1:t*9,:)=o2-o1;
        elas((t-1)*9+1:t*9,:)=omega((t-1)*9+1:t*9,:).*(1./shares*mp');
        markup((t-1)*9+1:t*9,1)=-inv(omega((t-1)*9+1:t*9,:).*(own_dum*own_dum'))*shares;
end
ivmelas(:,:,i)=median(reshape(elas',9,9,47),3)';
mc_hat=mean_price-markup;
company_post=company.*(1-(company==4))+(company==4)*3;
post_dum=cr_dum(company_post);

 options=optimset('Tolfun',1e-5,'TolX',1e-10,'Display','on');
flag=zeros(1,47);

first=1;
for j=1:47
	last=j*9;
    expall=exp(meany(:,j)-alpha*mean_price(first:last));
	[p_star(first:last,1),y,x]=fsolve(@logitbert_eq,mean_price(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
    flag(:,j)=x==1;
    first=last+1;
end
ivflag(i,:)=flag;
c=reshape(100*((p_star-mean_price)./mean_price),9,47)';
check=(ivflag(i,:)'+(sum(reshape(mean_share,9,47))<1)'==2);
ivpchange(i,:)=median(c(check,:));
end





%Test
%This evaluated first order conditions at post merger averages.  This is
%done by setting the tolerance on fsolve to be extremely high and then
%looking at the function evaluation at the "optimum".  

load 'postoutside4'
ols_evaluated=zeros(47,9);
mean_price=data(:,1);
mean_share=data(:,2);
meany=log(reshape(mean_share,9,47))-log(ones(9,47)-kron(ones(9,1),sum(reshape(mean_share,9,47))));
options=optimset('Tolfun',1e500,'TolX',1e-10,'Display','on');
first=1;
sims=1000;
olspchange=zeros(sims,9);
for z=1:sims
        bs_beta_IV=beta_IV+chol(vcov_IV_neweywest)'*randn(size(beta_IV,1),1);
        alpha=bs_beta_IV(1)
        markup_multi = zeros(9*47,1);
        elas = zeros(9*47,9);
        omega = zeros(9*47,9);

    


        o1=alpha*mean_share_atavg*mean_share_atavg';
        o2=diag(alpha*mean_share_atavg);
        omega=(o2-o1);
    
    for t=1:47
        shares=mean_share((t-1)*9+1:t*9);
     	mp=mean_price((t-1)*9+1:t*9);
        own_dum=pre_merger_dum((t-1)*9+1:t*9,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*9+1:t*9,:)=o2-o1;
        elas((t-1)*9+1:t*9,:)=omega((t-1)*9+1:t*9,:).*(1./shares*mp');
        markup((t-1)*9+1:t*9,1)=-inv(omega((t-1)*9+1:t*9,:).*(own_dum*own_dum'))*shares;
    end
    ivmelas(:,:,i)=median(reshape(elas',9,9,47),3)';
    mc_hat=mean_price-markup;
    company_post=company.*(1-(company==4))+(company==4)*3;
    post_dum=cr_dum(company_post);

    options=optimset('Tolfun',1e-5,'TolX',1e-10,'Display','on');
        
    for i=1:47
        last=i*9;
        expall=exp(meany(:,i)-alpha*post_mean_price(first:last));
        [p_star(first:last,1),y,ivflag(i,1)]=fsolve(@logitbert_eq,post_mean_price(first:last),options,expall,alpha,... 
        mc_hat(first:last), post_dum(first:last,:));
        first=last+1;
    end
    olspchange(i,.)=mean(p_star)';
end
